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We report a combined experimental and simulation study of deformation-induced diffusion in 
compacted two-dimensional amorphous granular pillars, in which thermal fluctuations play negligible 
role. The pillars, consisting of bidisperse cylindrical acetal plastic particles standing upright on a 
substrate, are deformed uniaxially and quasistatically by a rigid bar moving at a constant speed. 

The plastic flow and particle rearrangements in the pillars are characterized by computing the best- 
ht affine transformation strain and non-affine displacement associated with each particle between 
two stages of deformation. The non-affine displacement exhibits exponential crossover from ballistic 
to diffusive behavior with respect to the cumulative deviatoric strain, indicating that in athermal 
granular packings, the cumulative deviatoric strain plays the role of time in thermal systems and 
drives effective particle diffusion. We further study the size-dependent deformation of the granular 
pillars by simulation, and find that different-sized pillars follow self-similar shape evolution during 
deformation. In addition, the yield stress of the pillars increases linearly with pillar size. Formation 
of transient shear lines in the pillars during deformation becomes more evident as pillar size increases. 

The width of these elementary shear bands is about twice the diameter of a particle, and does not 
vary with pillar size. 

PACS numbers: 45.70.-n, 47.57.Gc, 83.50.-v, 83.80.Fg 


I. INTRODUCTION 

Disordered materials such as metallic glasses can 
exhibit highly localized deformation and shear band 
formation[l, 2]. Most experiments on these systems, 
however, use loading geometries in which there are free 
boundaries and inhomogeneous strains, while simulations 
have typically focused on systems with periodic boundary 
conditions under homogeneously-applied shear strain. To 
understand at a microscopic level the effects of loading 
geometry on the macroscopic mechanical response, it is 
useful to study a disordered system in which individual 
particles can be imaged and tracked as they rearrange 
under an applied load. Here we introduce a granular 
packing-a packing of discrete macroscopic particles for 
which thermal agitation plays a negligible role [3, 4]-in 
a pillar geometry commonly used for mechanical testing 
of metallic glasses. We combine experiment and simula¬ 
tion to study the response of these two-dimensional (2D) 
pillars to athermal, quasistatic, uniaxial compression. 

One question of interest is how the mechanical response 
of the pillar depends on pillar size. We find that the pillar 
shape evolves under load in a self-similar fashion, so that 
the shape of the pillar at a given strain is independent 
of system size. We also find that as the pillars deform, 
the strain rate localizes into transient lines of slip, whose 
thickness of a few particle diameters is independent of 
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system size. Thus, the system is self-similar in shape at 
the macroscopic scale, but, surprisingly, its yielding is 
not self-similar at the microscopic scale. 

A second question concerns the random motions of 
particles as they rearrange under inhomogeneous load¬ 
ing conditions. Because particles jostle each other, they 
display diffusive behavior in homogeneously sheared sys¬ 
tems that are devoid of random thermal fluctuations [5] . 
Recently, crystal nucleation and growth were observed 
in situ in mechanically fatigued metallic glasses at low 
temperature [6]. Crystallization is typically thought 
to require diffusion. Therefore, it was suggested that 
the “shear transformation zones” (STZs) [2] should 
be generalized to “shear diffusion transformation zones” 
(SDTZs) [6], to reflect the contributions of random mo¬ 
tions driven by loading, even under inhomogeneous con¬ 
ditions. Our amorphous granular pillar is an athermal 
system as far as the macroscopic particles are concerned 
(effective vibrational temperature « 0), so our experi¬ 
ment and simulations can examine how inhomogeneous 
loading affects particle motion. We find that the idea of 
load-induced diffusion can be generalized to inhomoge¬ 
neous loading by replacing time with the cumulative de¬ 
viatoric strain, and the mean-squared displacement with 
the mean-squared displacement of a particle relative to 
the best-fit affine displacement of its neighborhood {i.e. 
the mean-squared non-affine displacement [7] ). With this 
generalization, we observe that the mean-squared non- 
affine particle displacement crosses over from ballistic to 
diffusive behavior as a function of the cumulative devia¬ 
toric strain. 
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The article is organized as follows. In section II, we 
describe the experimental and simulation setup, as well 
as the simulation methodology, of 2D amorphous gran¬ 
ular pillars under uniaxial and quasistatic deformation. 
Section III describes the results of our combined experi¬ 
ments and simulations on the deformation of a 2D gran¬ 
ular pillar containing 1000 particles. In section IV, we 
discuss the exponential crossover of non-affine particle 
displacement from ballistic motion to diffusion with re¬ 
spect to cumulative deviatoric strain. Section V presents 
our simulation results on the size-dependent deformation 
of large 2D granular pillars. Then we conclude the article 
in section VI. 


II. METHODS 

The compacted 2D amorphous granular pillars in our 
study consist of 50-50 mixture of bidisperse cylindrical 
particles (grains) standing upright on a substrate. A top- 
view of the schematic setup is shown in Fig. 1. The pil¬ 
lars have aspect ratio Hq/Wq « 2, where Hq and Wq are 
the original height and width of the pillars respectively. 
In our experiment, the cylindrical granular particles are 
made of acetal plastic. The diameter of the large grains 
in the pillars, denoted by D, is 1/4 inch (0.635 cm), while 
for the small grains the diameter d has the value of 3/16 
inch (0.47625 cm). The ratio of diameter between large 
and small grains is therefore D/d = 4/3. Both types 
of grains are 3/4 inch (1.905 cm) tall. The masses for 
the large and small grains are 0.80 gram and 0.45 gram 
respectively. The pillars are confined between a pair of 
parallel bars. The bottom bar is static while the top bar 
deforms the pillars uniaxially with a slow, constant speed 
Vc = 1/300 inch per second (0.0084667 cm/sec). The 
force sensors connected to the bars measure the forces 
on the top and the bottom bars, and the trajectory of 
each particle in a pillar is tracked by a high-speed cam¬ 
era mounted above the pillar. The basic parameters in 
our simulation, including the size and mass of the grains, 
as well as the velocity of the bars, are the same as in 
the experiment. Further experimental details will be de¬ 
scribed in an upcoming paper [8]. 


A. Packing Generation Protocol 

Properly prepared initial configurations are crucial for 
the study of the mechanical properties of amorphous 
solids. In our experiment, 50-50 random mixture of bidis¬ 
perse grains are compacted to form a pillar with aspect 
ratio 2 to 1. To facilitate direct comparison between ex¬ 
periment and simulation, for small-sized pillars (num¬ 
ber of grains in the pillar N = 1000), the simulation 
initial conditions are taken from the experimental data, 
which was then relaxed in simulation to avoid particle 
overlapping resulted from measurement error. For large¬ 
sized pillars, which can only be studied by simulation, we 


compression speed Vc 



-► static bar 

Wo 


FIG. 1. Top-view of the experimental/simulation setup. A 
two-dimensional pillar of granular particles on a frictional sub¬ 
strate are deformed quasistatically and uniaxially by a rigid 
bar from one side. The direction of gravity is perpendicular 
to the substrate. The compacted, disordered granular pack¬ 
ing consists of 50-50 mixture of bidisperse cylindrical-shape 
grains. The ratio of radius between large and small grains 
is 4:3. The aspect ratio of the pillar, defined as the initial 
height of the pillar {Ho) divided by the initial width {Wo), is 
2:1. The pillar is confined between two rigid bars placed at 
the top and bottom end of the pillar respectively. The top bar 
deforms the pillar with a constant speed Vc while the bottom 
bar is kept static. 


generate compacted, amorphous granular pillars through 
computer simulation, using the protocol described below. 
The particle area density in the simulation-generated pil¬ 
lar is controlled to be at the onset of jamming transition 
[9]. To generate the initial conditions, we assign the fol¬ 
lowing truncated Lennard-Jones potential with purely re¬ 
pulsive interaction to the large (L) and small (S) grains 


Ua0{r) 


e {(Tap/rf^ - 2 {(Jap/rf 


for r < Gap, 
for r> Gap, 
( 1 ) 


where the subscripts a, /3 denote L or S. The zero-force 
cut-off distances Gap are chosen to be the sum of radii of 
two particles in contact, namely ctll = D, ctls = 711/8, 
and tTss = 3D/4, where D is the diameter of a large 
grain. We note that this potential will only be used to 
generate the initial conditions of the granular packings, 
and is different from the particle interaction model we 
describe later for the deformation of the granular pillars. 

To create a disordered granular packing with 50-50 
mixture of N total number of large and small grains, 
a rectangular simulation box with dimensions A x 2A is 
initially created, where the width of the box A is chosen 
such that the initial particle area density, p = N/2h/, 
is slightly above the particle overlapping threshold. We 
then randomly assign the positions of the particle within 
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the simulation box, and subsequently use conjugate- 
gradient (CG) method to minimize the total potential 
energy of the system. Periodic boundary conditions are 
applied during this process. The particle positions are 
adjusted iteratively until the relative change of energy 
per particle between two successive CG steps is smaller 
than 10“^^. When this stage is reached, the pressure of 
the system is calculated using the following virial formula 


P = - 


1 


Vr- — 


i>j 


( 2 ) 


where A is the area of the simulation box, is the dis¬ 
tance between particles i and j. If the pressure is greater 
than zero, both dimensions of the simulation box will be 
enlarged by a fraction of 10“^, and the particles in the 
box will be mapped to the corresponding new positions 
in the enlarged box via affine transformation. CG en¬ 
ergy minimization will then be carried out on the new 
configuration. This iterative process stops when the cal¬ 
culated pressure of the system at the end of a CG run 
becomes smaller than The final configura¬ 

tion will be taken as the initial conditions of close-packed 
2D amorphous granular assembly. Vacuum space is then 
added on the lateral sides of simulation box to create a 
pillar with 2:1 aspect ratio. Calculation of radial distri¬ 
bution functions for different-sized pillars indicates that 
the structure of the amorphous assemblies generated fol¬ 
lowing the above procedures does not show noticeable 
size dependence. Comparison of the radial distribution 
functions computed for the experimental and simulation¬ 
generated initial conditions is shown in Fig. 2. 


B. Simulation Methodology 

We use the method of MD to simulate the quasistatic 
deformation of the 2D granular pillars. The simula¬ 
tion force model includes three components: the grain- 
grain interaction, the grain-bar interaction and the grain- 
substrate interaction. Each of these forces will be de¬ 
scribed in the following. 


1. Grain-Grain Interaction 

As illustrated in Fig. 3a, the interaction between 
two grains includes normal and tangential contact force, 
which are denoted by F„ and Ft respectively. Two grains 
experience a repulsive normal contact force if the dis¬ 
tance between the particle centers is smaller than the 
sum of their radii. For two smooth, elastic cylindrical 
particles with parallel axes, the normal contact force as 
determined by the Hertzian theory of contact mechanics 
is proportional to the indentation depth between the two 
particles [10]. For our granular particles, denote by 
and Tj the positions of particles i and j, and denote by 
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FIG. 2. Comparison of the radial distribution functions g{r) 
for experiment-derived and simulation-generated initial con¬ 
ditions computed using (a) small grains as the central parti¬ 
cles and (b) larger grains as the central particles are shown 
respectively. The distance r is scaled by the diameter D of 
the large particles. 



FIG. 3. (a) Illustrations of grain-grain interaction in the 

granular pillar. The contact force between two grains consists 
of normal repulsive contact force F„ and tangential shear con¬ 
tact force Ft. (b) Illustration of grain-substrate interaction. 
If the velocity of a grain i is non-zero, or the vector sum of 
the forces on the grain due to other grains and the bars is 
non-zero, the substrate will exert a frictional force f on the 
grain, the maximum value of which is migfj,, where mi is the 
mass of the particle, g is the gravity acceleration constant 
and g denotes the friction coefficient between the grain and 
the substrate. Likewise, if the angular velocity of the grain is 
non-zero or the torque on the grain due to other interactions 
is non-zero, the substrate will induce a frictional torque whose 
maximum magnitude is |T^^i| = ^miggRi, where Rt is the 
radius of the particle. 
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Tij = Vi — rj the distance vector between the two parti¬ 
cles, the indentation depth Sij is calculated as 

Sij — Ri “t“ Rj ; ( 3 ) 

where rij = |rij|. Ri and Rj are the particle radius of 
i and j respectively. 6ij will be zero if the two particles 
are not in contact. The normal contact force acting on 
the particle i by particle j is then given by 


^riij — t (4) 

where Uij = r^- /r^-, and kn is the normal contact stiff¬ 
ness. The corresponding normal contact force on particle 
j is given by Newton’s third law, namely, = —F„ . 
In Hertzian theory of contact mechanics [10], the con¬ 
stant kn for two cylinders in contact can be calculated 
as 


kn = (5) 

where I is the height of the cylinders. E* is the normal¬ 
ized contact elastic modulus, which is computed from the 
respective elastic modulus of the two cylinders, Ei and 
E 2 , and their Poisson’s ratios, and i' 2 - 


1 I 1 - 

El E2 


( 6 ) 


The existence of a friction force between two particles 
in contact is a characteristic feature of granular materi¬ 
als. Appropriate modeling of contact friction is crucial to 
the study of granular dynamics. The tangential frictional 
force between two grains in contact can be very compli¬ 
cated in reality [11]. We adopt the history-dependent 
shear contact model initially developed by Cundall and 
Strack [12]. This well-tested model has been used by 
many others to model the dynamics of granular assem¬ 
blies [11, 13-21]. The essence of this model is to keep 
track of the elastic shear displacement of two particles 
throughout the lifetime of their contact, and applying 
the Coulomb elastic yield criterion when the displace¬ 
ment reaches a critical value. Our implementation of the 
Cundall-Strack model follows Silbert et al. [11]. Specifi¬ 
cally, the tangential contact force between particle i and 
j is calculated as: 


Ft.^- =(7) 

where the shear displacement is obtained by inte¬ 
grating the tangential relative velocities of the two par¬ 
ticles during the lifetime of their contact [11]. Here kt is 
the tangential contact elastic modulus. It is taken to be 
proportional to the normal contact stiffness kn- Follow¬ 
ing Silbert et al., we choose kt = ^kn- Previous studies 
showed that the dynamics of system is relatively insen¬ 
sitive to this parameter [11], which is confirmed by our 
own simulation. 

To model the elastic yield of shear contact, the mag¬ 
nitude of Ut-^ is truncated to satisfy the Coulomb yield 


criterion [Ft^.j < j/r^F^. j, where is the friction coef¬ 
ficient between the grains. 

The tangential contact force will induce torques on the 
two grains in contact, as given by 

Tij = ~2^ij ^ (^) 

Here T^- is the torque exerted by grain j on grain i due 
to the tangential contact force Fj^^.. 


2. Grain-Bar Interaction 

The grain-bar interaction is modeled in a similar way 
to the grain-grain interaction. The bar is essentially 
treated as a rigid grain with infinitely large radius. When 
a grain comes in contact with a bar, the grain can experi¬ 
ence normal and shear contact force induced by the bar, 
and the shear contact force is also calculated by tracking 
the elastic shear displacement between the grain and the 
bar. The motion of the moving bar is not affected by the 
grains. The static bar at the bottom side of the pillar is 
always static, while the top bar deforms the pillar at a 
constant speed Vc- Compared to grain-grain interaction, 
the interaction parameters between the grains and the 
bar is slightly modified. Since the bars are modeled as 
rigid bodies that cannot be elastically deformed, it means 
that the elastic modulus of the bars is considered to be 
infinite. Consequently, the effective interaction modu¬ 
lus E* between the bars and the grains, based on Eq. 6, 
is twice as large as that between the grains. Therefore, 
from Eq. 5, the normal interaction stiffness between the 
bars and the grains is twice as large as that between the 
grains, i.e., (grain-bar) = 2fc„(grain-grain). Since the 
shear modulus of contact kt is proportional to kn, we 
have fct (grain-bar) = 2fct(grain-grain) as well. 


3. Grain-Substrate Interaction 

The effect of the substrate on the grains is determined 
after all the forces and torques on each grain due to other 
grains and bars have been determined. The substrate can 
induce both frictional force and torque on the grains, as 
illustrated in Fig. 3b. If a grain is initially static, unless 
the magnitude of total force due to other grains/bars 
is larger than the maximum frictional force that can be 
exerted by the substrate [fi] = rmgfJ., the substrate fric¬ 
tional force will cancel out other forces on the grain and 
the particle will continue to have zero velocity. Here rrii 
is the mass of the grain i, g is the gravitational accelera¬ 
tion and g denotes the frictional coefficient between the 
grains and the substrate. In another case, if the veloc¬ 
ity of the grain is non-zero, the substrate will induce a 
frictional force opposite to the direction of particle mo¬ 
tion, with magnitude = rriigfi. A similar algorithm 
applies to the rotational motion of a particle. An initially 
static grain will not start to rotate unless the torque due 
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to other interactions surpasses the maximum substrate- 
induced frictional torque |T^J = ^rriig^iRi, where Ri is 
the radius of the cylindrical-shape particle. The prefac¬ 
tor I is based on the assumption that frictional force is 
evenly distributed on the circular contact interface be¬ 
tween a cylindrical-shape grain and the substrate. If 
the angular velocity of the grain is non-zero, a frictional 
torque 

2 

=--nuggRiLj^, (9) 

will slow down the rotational motion of the particles, 
where dJi = u}i/\u}i\ and uii denotes the angular velocity 
of particle i. 


two numbers, they cannot be damped further but instead 
oscillate around zero, which is clearly a numerical arti¬ 
fact. To work around this issue, we introduce two small 
parameters 

= gfj.6t, (11) 

Olii 

such that when |vi| < and | J2j ^ are 

both satisfied, the velocity and total force on the particle 
will be set to zero. Here F^- is the force of particle j on 
particle z, and is the force of the bars on particle 
i. Similarly, for the rotational motion, if \tjJi\ < and 
I Tij -I- T)’®'*'! < ^rriiggRi, the angular velocity and 
total torque of the particle is set to zero. 


4- Equations of Motion 


C. Choice of Simulation Model Parameters 


After all the forces and torques on an individual grain 
are determined, they are summed up and the velocities 
and angular velocities of the grains are then updated ac¬ 
cording to Newtonian equations of motion: 


rrii 


dt'^ 


= F,; 


, d(jJi ^ 

h—r^ = Ti 
dt 


( 10 ) 


where Fj and T* are the total force and torque on the 
particle i respectively. A = ^rtiiR^ is the moment of 
inertia for grain i. The standard velocity Verlet integra¬ 
tor is used to update the positions and velocities of the 
particles, while a finite difference method is used to inte¬ 
grate the first-order differential equation for the angular 
velocities. 

There is a subtle numerical issue that must be ad¬ 
dressed when modeling velocity and angular velocity 
changes of the particles in the presence of the damp¬ 
ing effects of a frictional substrate. In numerical inte¬ 
gration of equation of motion, time is discretized into 
small timesteps with each timestep being a small incre¬ 
ment 6t. To complete the simulation within a reasonable 
time frame, St cannot be too small, which means that 
the changes of velocity and angular velocity of the grains 
due to the substrate induced force and torque within a 
timestep are not infinitesimal. Hence, the motion of par¬ 
ticles might not be able to be brought to a halt by the 
substrate — the velocity and angular velocity of the parti¬ 
cles could oscillate around the zero. Consider, for exam¬ 
ple, a stand-alone cylindrical grain with initial velocity 
Vi and angular velocity Ui. Without other interactions, 
the substrate will induce friction jf^l = rriigg, and fric¬ 
tional torque = 2migg/3 on the grain, which slows 

down the translational and rotational motion of the grain 
respectively. According to the equations of motion in 
Eq. 10, the translational and rotational acceleration will 
be Qv = gg, and a^, — Agg/(ZRi), with Ri being the 
radius of particle i. Hence, within a timestep St, the 
change of velocity or angular velocity is a finite number: 
Sv = ggSt, Suj = 4ggSt/{3Ri). If the velocity or an¬ 
gular velocity have been damped to values below these 


The independent parameters in the interaction model 
of our simulation include the grain-grain stiffness 
grain-grain friction coefficient gg, grain-substrate fric¬ 
tion coefficient g, and the timestep for integration of 
equations of motion St. Among these parameters, g has 
been experimentally measured to be around 0.23. Hence 
g = 0.23 will be adopted in our simulations. The grain- 
grain friction coefficient gg is unknown. We have carried 
out simulation using multiple values of gg, and the re¬ 
sults indicate that choosing gg = 0.2 can achieve good 
match between the experiment and simulation. Due to 
the quasistatic nature of deformation by the moving bar 
on the pillars, the increment of force on a grain by the 
bar within one timestep St must be much smaller than 
the maximum static friction by the substrate on a grain, 
namely 


2knVcSt <C rriigg, (12) 

where Vc is the speed of the top moving bar. Hence, the 
smaller the value of St, the higher the value of that 
can be adopted in simulation. While there is no physi¬ 
cal reason for a lower bound of St, smaller St results in 
an increased time span to complete simulation. Realistic 
consideration leads to our choice of St = 10“® second. 
The upper bound of allowed calculated from Eq. 12 
is considered to be smaller than the real contact stiffness 
of two particles in experiment. For this reason, we have 
systematically studied the influence of on the simula¬ 
tion results in a small-sized pillar containing 1000 grains. 
The relatively small sized pillar allows us to use St = 10“® 
second and thus access a wider range of A:„, from /c„ = 
1 N/mm to kn = 100 N/mm. The results indicate that 
the statistical behaviors of deformation dynamics, such 
as flow stress and particle-level deformation characteris¬ 
tics, are not significantly influenced by the value of the 
kn- We therefore choose kn = 10 N/mm and St = 10“® 
in our simulation. 

The results of our study will be expressed in terms of 
several characteristic units. Length will be expressed in 
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the diameter of the large grains D or the radius R = D/2. 
The unit of velocity will be the bar speed Vc and the unit 
of time will be R/vc, which is the time it takes for top 
bar to move over a distance equal to R. The units for 
force and stress will be mg^, mgg/D respectively, where 
for convenience, we will use the symbol m to denote the 
mass of a large grain, mgg is thus the minimum force 
to induce the translational motion of a stand-alone large 
grain and mggi/D is the corresponding averaged stress of 
the bar on the grain. 


III. COMBINED EXPERIMENT AND 
SIMULATION ON DEFORMATION OF 
SMALL-SIZED PILLARS 

Deformation of an = 1000 pillar has been studied by 
both experiment and simulation. The experimental ini¬ 
tial particle arrangement in the pillar is the same as those 
depicted in Fig. 1. To facilitate comparison between ex¬ 
periment and simulation, our parallel simulation of pil¬ 
lar deformation uses the experimentally measured initial 
conditions, which were further relaxed in simulation to 
avoid particle overlapping resulted from measurement er¬ 
ror. When the pillar is deformed by the moving bar, the 
strain of deformation e is defined as the change of pillar 
height AH divided by the original height of the pillar 
Hq, namely, e = AH/Hq. The deformation stress a is 
calculated as the normal force on the top moving bar Fbar 
divided by the maximum width of the pillar near the top 
edge W, namely cr = Fh^r/W. 

Fig. 4 shows the experimental and simulation stress- 
strain curve of the N = 1000 pillar. The measured stress 
shows yielding behavior when the deformation strain ex¬ 
ceeds a very small value Sy. From our simulation, we 
find that the yield strain Sy in general becomes smaller 
as the grain-grain stiffness or the packing density of 
the pillar is increased. The yield stress ay however shows 
little dependence on The parameter that affects ay 
most was found to be the grain-grain friction coefficient 
g.g. In the range of g,g we have studied {g,g from 0 to 0.3), 
ay increases monotonically with the increase of fig. The 
simulation results presented in this paper use fig = 0.2, 
which was found to achieve overall good match between 
the experiment and simulation. 

In Fig. 4, we label several stress/strain values and cali¬ 
brate the corresponding particle-level structural changes 
in the pillar. The experimental and simulation results are 
then compared side-to-side in Fig. 5. Fig. 5a shows the 
mean particle velocity field in the pillar at six different 
stages of deformation. The mean velocity of a particle 
i, denoted by Vi{t,At), is calculated as the average dis¬ 
placement magnitude of the particle from current time t 
to a later time t + At, 


Vi{t,At) = \ri{t + At) - r^{t)\/At, (13) 

where the value of time interval At is chosen to be 2/15 
R/vc for the present purpose. Vi{t, At) contains informa- 



FIG. 4. Comparison between the (a) experimental and (b) 
simulation stress-strain curves for the deformation of a A = 
1000 granular pillar. The compressing stress is measured in 
units of mgfi/D, while the strain is computed as the change 
of pillar height (AH) divided by the original height of the 
pillar Ho- The special labels (1-6) indicate the stress strain 
values at which deformation characteristics in the pillar will 
be compared side to side between experiment and simulation. 


tion of the absolute amount of displacement of the parti¬ 
cle i within At. As shown in Fig. 5a, the mean velocities 
of the particles near the moving bar are close to v^, which 
is expected as the pillar is deformed quasistatically by the 
bar. The mean velocity of a particle in general becomes 
smaller as the particle is further away from the moving 
bar. At the early stages of deformation, particles at the 
bottom part of the pillar have not moved and therefore 
have zero values of u. A sharp boundary between the 
moving and non-moving regions of the pillar often forms 
along the the direction that is roughly 45 degree to the 
direction of uniaxial deformation. 

In the simulations we have access to detailed informa¬ 
tion on the inter-particle interactions. In Fig. 6 we plot 
the grain-grain normal force Fn, tangential force Ft and 
substrate-induced force frictional force / on the particles 
at six stages of deformation corresponding to the special 
labels in Fig. 4. Comparing Fig. 6a with Fig. 6b and 
Fig. 6c, we find that F^ is in general much larger than 
Ft, which is further larger than /, namely F^ ^ Ft ^ f . 
In particular, Fig. 6a shows that particles with large Fn 
are connected with force chains. The magnitude of forces 
in these force chains is higher for particles residing in the 
interior the pillar. This indicates that the stress in the 
pillar is rather inhomogeneous, with larger stresses in the 
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FIG. 5. Comparison between experiment and simnlation of the particle velocity v, deviatoric strain rate J 2 and non-afiine 
displacement during deformation of a = 1000 granular pillar. The six stages of deformation (1-6) correspond to the 

stress and strain values labeled in Fig. 4. Within each subplot (a), (b) and (c), the top panel corresponds to the experimental 
result, while the bottom panel corresponds to the simulation result, (a) Velocities of the particles in the pillar. The magnitude 
of the displacement of a particle from the current position after time interval At = {2/15)R/vc is divided by At to obtain 
the average velocity across the time interval, (b) Deviatoric strain rate J 2 for each particle. J 2 is computed by comparing 
the current configuration of a particle and its neighbors with the configuration after At, using neighbor sampling distance 
Rc = 1.25D. J 2 is measnred in the unit of VcjR- (c) Non-afBne displacement for each particle in the pillar. The 

procedures for calculating are discussed in the main text. 


interior region of the pillars than close to the surface. 

We further look at the rearrangement of particles in the 
pillar by defining a neighbor sampling distance i?c, and 
calculate the affine transformation strains and non-afhne 
displacements of the particles with respect to their neigh¬ 
bors within Rf.. The value of is chosen to be 1.25Z1, 
which roughly corresponds to the average first nearest- 


neighbor distance of the particles in the pillar, as can be 
seen from the computed radial distribution functions in 
Fig. 2. A particle j is considered to be the neighbor of 
a particle i if their distance is smaller than i?c, which is 
illustrated in Fig. 7. The configurations of the particle 
i and its neighbors at a given time t and a subsequent 
time t -I- At will then be used to compute the best-fit 
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FIG. 6. Forces in the granular pillar during deformation as obtained from simulation. The special labels (1-6) correspond 
to the deformation stages labeled in Fig. 4b. (a) Grain-grain normal force Fn\ (b) grain-grain tangential force Ft and (c) 
grain-substrate friction force /. The normal and tangential forces are measured in the unit of mgfj,, which is the largest possible 
value of substrate induced friction on a large grain. The substrate friction forces are measured in the unit of rriigg. 


local affine transformation matrix J and the non-affine 
displacement associated with particle i, using the 
method introduced by Falk and Langer [7, 22]. Specif¬ 
ically, i is obtained by calculating the best affine 

transformation matrix that minimizes the error of de¬ 
formation mapping: 

At) = ;^min + At) - , 

* ’ jeNi 

(14) 

where rji{t) = — is the distance vector between 

particles j and i at time t. rji{t -|- At) is the distance 
vector at a later time t -|- At. The summation is over the 
neighbors of particle i at time t, whose total number is 
given by Ni. The best-fit affine transformation matrix 
Ji(t, At) is usually non-symmetric due to the presence of 
rotational component. A symmetric Lagrangian strain 
matrix tj^ can be calculated from as 

= (15) 

where I is an identity matrix. The hydrostatic invariant 
is then computed from as 

rir = ^Trr,,. (16) 


The shear (deviatoric) invariant is then given by 

Vt = y^Tr(r7, (17) 

Hereafter we will refer to 77 ® (t, At) as the deviatoric strain 
associated with the particle i from t to t -|- At. The 
deviatoric strain rate, denoted by J 2 , is the normalization 
of 77 ® (t. At) with respect to At: 

J2{t, At) = 77®(t, At)/At (18) 

Fig. 5b-c shows the computed deviatoric strain rate J 2 
and for each particle in the pillar at six different 

stages of deformation, where the experimental and sim¬ 
ulation results are compared side to side. J 2 {t,At) and 
Il^;„(t, At) are computed using At = {2/15)R/vc, which 
is the same as the value of At used for computing the 
mean velocities of the particles in Fig. 5a. Comparing 
Fig. 5b with Fig. 5a, it can be seen that large values of 
deviatoric strain rate occur at places where the gradient 
of mean velocity, and hence the gradient of particle dis¬ 
placement, is large, which is understandable as strain is 
a measure of displacement gradient. One can also no¬ 
tice from Fig. 5b the presence of thin shear lines in the 
pillars, where particles with large deviatoric strain rate 
reside. The width of these shear lines is about twice the 
diameter of the particles. These shear lines largely cor¬ 
respond to the moving boundary between the deformed 
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FIG. 7. (a) Illustration of a particle (colored in green) 

and its neighbors (colored in black) within a cut-off distance 
Rc = 1.25Z) at an initial reference conhguration. (b) The 
same set of particles at a later stage of deformation. We seek 
to find the best-affine transformation matrix J that maps the 
coordinates of the particles illustrated in (a) to those in (b). 
This optimization procedure also gives the non-affine displace¬ 
ment Umin associated with the central (green) particle, and 
the deviatoric strain rj” in the neighborhood, as discussed in 
the main text. 


CM 




CM E 

Q 


CC 



10 “^ 10 "^ 1 10 
deviatoric strain r]‘ 


Simulation 


Experiment 


FIG. 8. Fitting of with respect to deviatoric strain 

rf for both (a) experiment and (b) simulation. The data are 
fitted to an exponential crossover equation from quadratic to 
linear scaling (see the main text for details). In (a), the green 
curve is the best-fitting curve for the experimental data. For 
(b), as the fitting results is nearly identical to (a), we plot the 
same fitting curve as in (a) to demonstrate the closeness in 
the fitting result. 


and undeformed regions in the pillar. The presence of 
such shear lines will appear clearer as pillar size increases, 
which will be discussed in the later part of the article. 


Comparing the profile in Fig. 5c with deviatoric 
strain rate Ji in Fig. 5b, it is clear that particles with 
larger values of correlated with larger values 

of J 2 , and hence also deviatoric strain rf (Eq. 18). The 
deviatoric strain 77 ® reflects the local shear component of 
affine deformation (shape change), while measures 

additional particle displacement with respect to its neigh¬ 
bors that cannot be described by mere shape change. The 
positive correlation between and rf is understand¬ 

able because the larger the value of 77 ® (which usually 
drives plastic deformation), the error of describing local 
particle rearrangement in terms of purely shape change, 
which is the definition of more likely to be 

larger. 


IV. LOCAL DEVIATORIC STRAIN DRIVEN 
PARTICLE DIFFUSION 

The positive interdependence between and 77 ® 

motivates us to map out their correlation quantitatively. 
Starting with an initial configuration of the pillar at time 
t that corresponds to deformation strain e = Vct/Ho, 
we fix the neighbor sampling distance Rc = 1.25D and 
calculate 77 ®(t, At) and At) for each particle in 

the pillar using a logarithmic series of time interval 
At e [2,4, 8 ,..., I28J/15 R/vc- This procedure is then 
repeated for at least eight values of initial times t equally 
spaced by “^R/vc- We then plot all the calculated val¬ 
ues of At) with respect to 77 ®(t, At) on a single 

plot, using logarithmic axes for both and 77 ®. The 

results of experiment and simulation are shown together 
in Fig. 8 . From Fig. 8 , it can be seen that while for 
a given specific value of 77 ®, the possible values of 
are scattered, the existence of statistical correlation be¬ 
tween and 77 ® is apparent. We find that in the range 
of small values of 77 ®, scales quadratically with 77 ®, 

which gradually transits to linear scaling at larger values 
of 77 ®. This is reminiscent of the scaling relationship be¬ 
tween the growth of mean squared displacement (MSD) 
for a thermally diffusive particle and time t, which is often 
explained pedagogically by an unbiased random walker. 
Indeed, we find that, by considering D'^-^ as MSD, and 
deviatoric strain 77 ® as time, the data in Fig. 8 can be fit 
very well using the following equation that describes the 
exponential crossover of a thermal particle from ballistic 
to diffusive motion, expected for a Langevin particle with 
no memory [23]: 

= 40??'' - 4077® [1 - exp(- 77 ®/ 77 ®)], (19) 

where on the left hand side of the above equation, the cal¬ 
culated is scaled by R? to render it dimensionless. 
0 is the dimensionless effective diffusivity while 77 ® takes 
meaning of “crossover deviatoric strain”. Our fitting of 
the data gives 0 = 0.3, 77 ® = 0.049 for the experiment, 
and 0 = 0.3, 77 ® = 0.05 for the simulation. 

The analogy between and MSD, and between 77 ® 
and time t, may have deep implications. describes 
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the mean-squared non-affine displacement of a particle 
with respect to its neighbors and can be naturally iden¬ 
tified as an analogy to MSD. The analogy between de- 
viatoric strain 77 ® and time t implies that, for the gran¬ 
ular packings, where there is no thermal agitation and 
the system is deformed heterogeneously, the cumulative 
deviatoric strain plays the role of time and drives effec¬ 
tive particle diffusion. Argon had originally used bub¬ 
ble raft deformation to illustrate the concept of shear 
transformation zone (STZ) [24, 25], which emphasizes the 
affine part of localized stress-driven processes. Recently, 
Wang et al. found that cyclic mechanical loading can in¬ 
duce the nano-crystallization of metallic glasses well be¬ 
low the glass transition temperature [ 6 ], resulting from 
stress-driven accumulation of non-afRne displacement of 
the atoms in the sample. The concept of shear diffu¬ 
sion transformation zones (SDTZs) was proposed by the 
authors to explain the experimental results and to em¬ 
phasize the diffusive character of STZs. Our results lend 
support to the concept of SDTZ by showing that, even in 
amorphous granular packings, where there is no thermal- 
driven diffusion at all, if the accumulated local deviatoric 
strain is large enough (above a few percent strain), the 
non-affine displacement of a particle with respect to its 
neighbors crosses over to the diffusive limit. This sug¬ 
gests that SDTZ may be a key concept for a broad range 
of amorphous solids. 

The analogy between local cumulative shear transfor¬ 
mation strain in athermal amorphous solids and time in 
thermal systems for particle diffusion may be rational¬ 
ized by a “space-time equivalence” argument, as follows. 
A finite temperature means temporally random mo¬ 
mentum fluctuations, for crystals and non-crystals alike. 
Even in crystals, such random momentum fluctuations 
(due to collision of multiple phonons) can drive the ran¬ 
dom walker behavior of a particle, if these temporal fluc¬ 
tuations can be significant compared to the potential en¬ 
ergy barrier. But in amorphous solids without sponta¬ 
neous temporal fluctuations, there will be nonetheless 
still another source of randomness, which is the local spa¬ 
tial structure and structural response of the amorphous 
solid. This is indeed what motivated the “heteroge¬ 
neously randomized STZ model” [26, 27]. In other words, 
even if two “Eshelby inclusions” at different locations 
of an amorphous solid transform by exactly the same 
transformation strain 77 , one reasonably would still ex¬ 
pect drastically different internal particles arrangements 
and rearrangements inside these zones. This ultimately 
is because the local strain 77 is just a coarse-graining 
variable, that represents a key aspect of the structural 
transformations of a kinetically frozen random cluster, 
but not all of its structural information. (This may not 
be true in simple crystals, where 77 may entirely capture 
the entire structure.) Such structural mutations beyond 
transformation strain are reflected in The fact 

that will accumulate linearly with strain at steady 

state means the structural mutations from generation to 
generation [26, 27] are largely non-repeating and essen¬ 


tially unpredictable, if starting from a spatially random 
configuration at the beginning, even when the stress con¬ 
dition driving these transformations remains largely the 
same. Our experiment and simulation on compressing 
amorphous granular pillars can thus be seen as a “spatial 
random number generator” with the initial configuration 
as the “random number seed”, in contrast to more well- 
known “temporal random number generator” algorithms; 
but both types of algorithms tend to give long-term un¬ 
correlated increments for the random walker. 


V. SIMULATION OF SIZE-DEPENDENT 
PILLAR DEFORMATION 

Having achieved good agreement between experiment 
and simulation for the N = 1000 pillar, we now take 
advantage of the fact that our simulation can treat 
mnch larger systems than experiment, to stndy the size- 
dependent deformation behavior of the granular pillars 
by simulation. Three large-sized pillars, denoted by 
N = 4000, N = 16000 and N = 64000, are deformed 
by the top bar moving at the same deformation speed Vc- 
The aspect ratio of the pillars (2 to 1) is fixed to be the 
same value of the N = 1000 pillar. As the initial packing 
density of the particles in the pillar is also the same, the 
initial width of the pillars Wq scales as '/N. 

We find the macroscopic shape evolution of the 
different-sized pillars is self-similar during deformation. 
At the same valnes of deformation strain e = AH/Hq, 
we extract the boundaries of the pillars, rescale them by 
the respective initial pillar width Wq, and plot them to¬ 
gether in Fig. 9. The rescaled boundaries of the pillars are 
nearly identical to each other. This also implies that, the 
width of the top edge of a pillars W divided by its origi¬ 
nal width Wg, is to a good approximation only a function 
of strain e but not the pillar size, namely W/Wq = x(e), 
where the scaling function y does not depend on the pil¬ 
lar width Wq. 

We also find that, the average flow stress of the pil¬ 
lars increases linearly with the initial pillar width Wq, 
as shown in Fig. lO(a-b). Mathematically, this can be 
expressed as (cr) oc Wq, where we define (cr) to be the av¬ 
erage flow stress for strain e between 0.05 and 0.2. This 
scaling behavior for the flow stress indicates that, for the 
2D disordered granular pillars, the behavior of “smaller 
is weaker” is exhibited. This is quite different from the 
deformation of free-standing metallic glass pillars, where 
“smaller is stronger” is the general trend [28, 29]. 

To understand the surprising size dependence of flow 
stress, we first look at the stress distribution in the pil¬ 
lars. In Fig. 5 we have shown that the grains in the in¬ 
terior region of the pillar experience larger inter-particle 
contact forces, resulting in larger local stress in the in¬ 
terior region of the pillar. The rate of increase for local 
stress as a function of distance to the lateral edges of the 
pillars is found to be very close for different-sized pillars. 
Such stress non-uniformity should also be reflected in the 


11 


a 


— N=4000 
— N=16000 

— N=64000 


b 




d e f 





FIG. 9. Self-similar evolution of pillar shapes during deformation of different-sized pillars. The boundaries of three pillars 
{N = 4000, N — 16000 and N — 64000) are rescaled and plotted together at the same strain value. 


local contact pressure between the moving bar and the 
pillar. Indeed, we find that the contact pressure is also 
spatially rather non-uniform. Fig. 10c shows that, the lo¬ 
cal contact pressure increases almost linearly from near 
zero at the edge of pillar to saturated values around the 
center of contact interface. The maximum values of local 
contact pressure scale roughly linearly with pillar width, 
consistent with the linear scaling of pillar flow stress. 

Since the pillars are deformed quasistatically, most of 
the deformation work on the pillars will be dissipated 
during plastic flow. The flow stress is therefore closely 
related to the dissipation of energy in the systems. We 
hence study how the energy dissipation in the pillars 
changes with pillar size. As the granular particles in the 
pillars stand on a substrate, two major mechanisms of 
energy dissipation during plastic flow can be identified: 
one is due to the grain-substrate friction and the other 
can be attributed to the grain-grain friction. The total 
external power input by the moving bar into the pillar, 
denoted by Pin, can be calculated as 

Pin = FhaiVc = a-Wvc. (20) 

We have shown that, compared at the same deformation 
strain e, both the flow stress a and pillar width W are 
proportional to the initial pillar width Wq. Hence, the in¬ 
put power by the external force scales quadratically with 
Wo, namely Pin oc Wg. As most of the input power will 
be dissipated in the plastic flow regime, the dissipated 
power should also scale with Wg . To study how the dissi¬ 
pated power is distributed between the substrate-induced 
friction and grain-grain friction, we compute the frac¬ 


tion of input power dissipated by the grain-substrate fric¬ 
tional force and study its size dependence. The amount 
of power dissipated by the grain-substrate friction force, 
denoted by Pg.g, can be calculated as 

Pg-s = rriigfiv,, (21) 

i 

where the particle mean velocity Vi has the same defini¬ 
tion as in Eq. 13, namely the average displacement of the 
particle i within a small time interval At. The fraction 
of power dissipated by the substrate-induced friction, de¬ 
noted by K, is then given by /c = Pg-s/Pin- We calculate 
the values of k for different sized pillars and plot them 
as a function of deformation strain in Fig. lOd. The re¬ 
sult indicates that k is statistically independent of pillar 
size. This allows us to conclude that the amount of input 
power dissipated by grain-substrate friction, Pg_s = nPin, 
also scales quadratically with pillar size Wg, and hence 
scales linearly with the number of particles in the pillar 
N. This effectively means that the number of particles 
participating in the plastic flow scales linearly with the 
total number of particles in the pillar, which is consistent 
with the self-similar shape evolution of the pillars. 

The calculated values of k in Fig. lOd indicate that 
the majority of deformation work is dissipated by the 
friction between the particles in the pillar and the sub¬ 
strate. Substrate friction therefore must play an impor¬ 
tant role in the linear increase of flow stress with re¬ 
spect to pillar width and the self-similar evolution of 
pillar shape, which have been shown to be consistent 
with each other. The granular pillars in our study are 
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FIG. 10. Size-dependent flow stress and dissipation of input power, (a) Stress strain curves for different sized pillars, (b) 
Linear scaling of average flow stress with respect to pillar width Wo. The average flow stress is computed for the range of strain 
between 0.05 and 0.2. (c) Local contact pressure p between the moving bar and the pillars as a function of position x along the 
contact interface, computed for different sized pillars at a representative value of deformation strain in the plastic flow regime. 
The position x is scaled by the width of the pillar W at the contact interface, (d) Fraction of input power dissipated by the 
grain-substrate friction as a function of deformation strain for different-sized pillars. 


not truly two-dimensional due to the presence of grain- 
substrate friction. This setup is however necessary for 
stable plastic flow of the uniaxially deformed granular 
pillars without cohesive interparticle interaction. With¬ 
out the grain-substrate friction, the deformation behavior 
of the granular pillars are expected to be quite different, 
and the size-dependent deformation behavior observed in 
this study {i.e. “smaller is weaker”) may no longer hold. 

If the macroscopic shape evolution of the pillars in our 
systems is self-similar, then how does the local yield¬ 
ing behavior vary with pillar size? We characterize the 
deformation-induced local structural change of the pil¬ 
lar by computing the deviatoric strain rate J 2 associated 
with each particle between two stages of deformation, us¬ 
ing the same methodology described earlier in the article. 
We find that, within a small amount of pillar strain, par¬ 
ticles with large values of J 2 organize into thin shear lines, 
which becomes more evident as pillar size increases, as 


shown in Fig. 11. These shear lines orient along the di¬ 
rection about 45 degree to the direction of uniaxial com¬ 
pression. Clearly, such shear lines form along the direc¬ 
tion of maximum shear stress. The sharpest shear lines 
appear predominantly at the moving boundary between 
the deformed and undeformed region in the pillars, as 
we have mentioned earlier when discussing the combined 
experimental and simulation study of small-sized pillars. 
A close-up view of these shear lines in Fig. 11 indicates 
that the width of the shear lines does not change as pil¬ 
lar size increases, maintaining a value about twice the 
diameter of a grain. We emphasize that these shear lines 
are transient in time. As deformation goes on, new shear 
lines will form elsewhere in the pillar, while the particles 
in the shear lines formed earlier may not accumulate sig¬ 
nificant amount of shear strain continuously. Evidence 
of such transient shear bands in granular materials was 
previously reported in the discrete element simulations 
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A/= 1000 A/= 4000 A/= 16000 A/= 64000 

FIG. 11. Deviatoric strain rate J 2 associated with each particle and zoom-in views of the transient shear lines in different-sized 
pillars. Four different-sized pillars are compared with each other, which contain 1000, 4000, 16000, 64000 grains respectively 
(from left to right). The regions in the pillars for zoom-in views are indicated by squares. For each pillar, the two configurations 
of pillars used for J 2 calculation are separated by time difference At = 8/15 R/vc- 




by Aharonov and Sparts [30] and Kuhn [31, 32]. Mal¬ 
oney and Lemaitre [33], and Tanguy et al. [34] observed 
transient lines of slip in their athermal, quasistatic simu¬ 
lation of 2D glasses of frictionless particles, and explained 
their formation in terms of elastic coupling and cascad¬ 
ing of shear transformation zones. The results of our 
combined experiment (Fig. 5) and simulation of uniax¬ 
ial, quasistatic deformation of 2D granular pillars clearly 
demonstrate the existence of such transient shear lines, 
which carry localized deformation in the granular pillars. 

The size-independent width of the transient shear lines 
is surprising since the overall macroscopic shape of the 
pillar is self-similar in systems of different sizes. Despite 
the self-similarity at the macroscopic scale, the system is 
not self-similar in how it yields at the microscopic scale. 
Since the slip lines are independent of system size, there 
must be more of them in larger systems, which is indeed 
observed in our simulation. Why the system chooses to 
be self-similar at the macroscopic scale but not at the 
microscopic scale is an interesting point for future study. 


VI. CONCLUDING REMARKS 

We have carried out combined experiments and sim¬ 
ulations of the quasistatic, uniaxial deformation of 2D 
amorphous granular pillars on a substrate. The simu¬ 
lation model developed in this article achieves excellent 
quantitative match to the experiment. In particular, we 
find that, for the granular packings, the non-affine dis¬ 
placements of the particles exhibit exponential crossover 
from ballistic motion to diffusion-like growth behavior 
with respect to local deviatoric strains. This result is a 
generalization to inhomogeneous loading of earlier obser¬ 


vations of stress-driven diffusion of particles in simulated 
2D molecular glasses under simple shear or pure shear 
in the thermal, quasistatic limit [5, 34-38]. Because in 
our study the “time” variable for diffusion, the best-fit 
deviatoric strain in a neighborhood, is a local measure 
of deformation and shear transformation, we expect that 
the non-afhne displacement should cross over from bal¬ 
listic to diffusive behavior in amorphous solids under any 
loading conditions. 

In metallic glass pillars, the apparent strength of the 
pillar and strain localization behavior depends on pillar 
diameter, manifesting so-called “size-dependent plastic¬ 
ity” behavior [27]. Often, “smaller is stronger” holds for 
metallic glasses [28, 29]. We have shown that for 2D 
granular pillars on a substrate, the frictional interaction 
between the granular particles and the substrate leads to 
the opposite size-dependent response, namely “smaller is 
weaker”. 

Finally, our combined experiment and simulation 
study clearly demonstrate that transient lines of slip form 
in quasistatically deformed amorphous granular pillars 
under uniaxial loading condition. These system-spanning 
shear lines carry localized shear transformations in 2D 
granular pillars, and their width shows no size depen¬ 
dence. Altogether, these results could have important 
implications for the plasticity and internal structural evo¬ 
lution of amorphous solids. 
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